# -*- coding: utf-8 -*-
"""
Created on Mon Apr 22 15:34:55 2019

@author: Soly Liang
"""
from numpy import *
def LogisticRegression(X,Y,coe=1E-5):
    # 牛顿法对率回归
    # 输入X和Y，维度分别为m*n和m*1
    # coe为终止条件，当梯度小于coe时停止迭代
    # 返回拟合到的参数theta,其中最后一个参数theta[-1]为b
    m,n=shape(X)
    X1=c_[X,ones(m)]    #用c_函数实现数组的横向拼接
    theta=zeros(n+1)    #在matlab中这样生成的是 (n+1)*(n+1)的矩阵，
                        #但是在numpy中为1*(n+1)的向量
                        #注意这里向量不分行向量和列向量的
    dL=ones(n+1)        #梯度，初始任意设置为较大向量
    while linalg.norm(dL)>coe:
        p1=1/(1+exp(-dot(X1,theta)))          #计算p1,对应于(3.23)式，也就是logistic函数
        dL=-dot(X1.T,Y-p1)                    #梯度,结果为(n+1)维向量
        H=dot(dot(X1.T,diag(p1*(1-p1))),X1)   #hessian矩阵，结果为(n+1)*(n+1)
        theta-=dot(linalg.inv(H),dL)          #更新theta
    return theta                              #返回最佳参数

def LogisticPred(x_test,y_test,theta):
    #对率回归预测并给出错误率
    m,n=shape(x_test)
    x=c_[x_test,ones(m)]
    y_pre=zeros(m)
    y_pre[dot(x,theta)>0]=1
    return sum(y_pre!=y_test)/m,y_pre
    

def load_iris(path):
    #载入鸢尾属(Iris)植物数据,特征数4，分类数3，样本数50+50+50
    #输入path-文件路径，字符串类型，注意路径中斜杠的地方要用双斜杠
    feature_name=['sepal_length','sepal_width','petal_length','petal_width']
    class_name,x,y=[],[],[]
    f=open(path,'r')
    for line in f:
        ls=line.strip('\n').split(',')
        if len(ls)<len(feature_name)+1:  #最后一行为空，剔除该情况
            continue
        x.append([float(i) for i in ls[:-1]])
        if ls[-1] not in class_name:
            class_name.append(ls[-1])
        y.append(class_name.index(ls[-1])+1)
    f.close()
    return array(x),array(y),feature_name,class_name

def load_wine(path):
    #载入wine数据，特征数13，分类数3，样本数59+71+48
    #输入path-文件路径，字符串类型，注意路径中斜杠的地方要用双斜杠
    x,y=[],[]
    f=open(path,'r')
    for line in f:
        ls=line.strip('\n').split(',')
        if len(ls)<1:  #最后一行为空，剔除该情况
            continue
        x.append([float(i) for i in ls[1:]]) #类号放在第一列的
        y.append(float(ls[0]))               #wine数据集中类号直接为'1','2','3'
    f.close()
    return array(x),array(y)

def Class2(x,y,n1,n0):
    #仅选择原数据中的两个类作为二分类问题
    #n1对应的类y=1，n0对应的类y=0
    m1=len(y[y==n1])
    m0=len(y[y==n0])
    x1=x[y==n1,:]
    x0=x[y==n0,:]
    return r_[x1,x0],r_[ones(m1),zeros(m0)]    

def CrossValid(x,y,k,ceo=1E-5):
    #k折交叉验证法,k=m时即为留一法
    if k>len(y):
        print('错误！划分数大于样本数，无法划分。')
        return []
    #按类别拆分数据
    y1=y[y==1]; x1=x[y==1,:]; m1=len(y1)
    y0=y[y==0]; x0=x[y==0,:]; m0=len(y0)
    #将原数据打乱
    random.shuffle(x1)
    random.shuffle(x0)
    #num1，num0分别为划分的各个子集中正例和反例样本数
    num1=(m1//k)*ones(k,dtype=int); num1[:m1%k]+=1
    num0=(m0//k)*ones(k,dtype=int); num0[:m0%k]+=1
    #各子集中正例和反例的索引号
    index1=k*[[]]
    index0=k*[[]]
    #在列表index1,index0中填入索引号
    p1=0; p0=0
    for i in range(k):
        index1[i]=list(range(p1,p1+num1[i]))
        p1=p1+num1[i]
        index0[i]=list(range(p0,p0+num0[i]))
        p0=p0+num0[i]
    #将num0和index0倒序
    num0=num0[::-1]
    index0=index0[::-1]
    #k份子集中，每一份分别作为测试集，其余作为训练集
    erro=[]
    for i in range(k):
        x_test=r_[x1[index1[i],:],x0[index0[i],:]]
        y_test=r_[ones(num1[i]),zeros(num0[i])]
        #其余k-1份子集的序号
        other1=list(set(range(m1))-set(index1[i]))
        other0=list(set(range(m0))-set(index0[i]))
        x_train=r_[x1[other1,:],x0[other0,:]]
        y_train=r_[ones(len(other1)),zeros(len(other0))]
        #训练
        theta=LogisticRegression(x_train,y_train,ceo)
        #测试
        erro1,_=LogisticPred(x_test,y_test,theta)
        erro.append(erro1)
    return mean(erro)
    

#主程序

#在鸢尾花数据集上进行训练和测试
path='iris.data'
x,y,feature_name,class_name=load_iris(path) #载入数据
x1,y1=Class2(x,y,1,3)      #在1,2,3中随意选择两个类别进行二分类计算
erro_10fold=[]
for i in range(10):
    erro_10fold.append(CrossValid(x1,y1,10))   #10折交叉验证法
erro_LOO=CrossValid(x1,y1,len(y1)) #留一法，k=m的特殊情况
print('iris数据集-10折交叉验证法估计的错误率为：')
print('平均：{:.2%}'.format(mean(erro_10fold)))
print('最小：{:.2%}'.format(min(erro_10fold)))
print('最大：{:.2%}'.format(max(erro_10fold)))
print('iris数据集-留一法估计的错误率为：{:.2%}'.format(erro_LOO))

#在wine数据集上进行训练和测试
path='wine.data'
x,y=load_wine(path) #载入数据
x1,y1=Class2(x,y,2,3)      #在1,2,3中随意选择两个类别进行二分类计算
erro_10fold=[]
for i in range(10):
    erro_10fold.append(CrossValid(x1,y1,10))   #10折交叉验证法
erro_LOO=CrossValid(x1,y1,len(y1)) #留一法，k=m的特殊情况
print('wine数据集-10折交叉验证法估计的错误率为：')
print('平均：{:.2%}'.format(mean(erro_10fold)))
print('最小：{:.2%}'.format(min(erro_10fold)))
print('最大：{:.2%}'.format(max(erro_10fold)))
print('wine数据集-留一法估计的错误率为：{:.2%}'.format(erro_LOO))